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Abstract 

We investigate, analytically and numerically, families of bright solitons in a system 
of two linearly coupled nonlinear Schrodinger/Gross-Pitaevskii equations, describ- 
ing two Bose-Einstein condensates trapped in an asymmetric double-well potential, 
in particular, when the scattering lengths in the condensates have arbitrary mag- 
nitudes and opposite signs. The solitons are found to exist everywhere where they 
are permitted by the dispersion law. Using the Vakhitov-Kolokolov criterion and 
numerical methods, we show that, except for small regions in the parameter space, 
the solitons are stable to small perturbations. Some of them feature self-trapping of 
almost all the atoms in the condensate with no atomic interaction or weak repulsion 
coupled to the self-attractive condensate. An unusual bifurcation is found, when the 
soliton bifurcates from the zero solution without a visible jump in the shape, but 
with a jump in the number of trapped atoms. By means of numerical simulations, 
it is found that, depending on values of the parameters and the initial perturbation, 
unstable solitons either give rise to breathers or completely break down into inco- 
herent waves ("radiation"). A version of the model with the self-attraction in both 
components, which applies to the description of dual-core fibers in nonlinear optics, 
is considered too, and new results are obtained for this much studied system. 
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1 Introduction 



Experimental observation of Bose-Einstein condensates (BECs) in trapped dilute gases 
[1,2,3,4] had opened new exciting possibilities for manifestations of nonlinear phenomena 
in various geometries. Indeed, in the mean-field approximation (which usually applies to 
a great accuracy), the order parameter, which can be identified with the single-atom wave 
function, obeys the nonlinear Schrodinger (NLS) equation with an external potential, also 
known as the Gross-Pitaevskii (GP) equation [5] . By varying the trap potential, the shape 
of the condensate can be tailored to sundry geometries, in particular, the shape of the 
condensate may be approximately circular, i.e., the condensate itself may be effectively 
two-dimensional (2D) or one-dimensional (ID, alias cigar-shaped) [6]. 

It was recently shown that not only the geometry of the condensate, but also the magni- 
tude and the sign of the scattering length, which determines interactions between atoms 
in the condensate and, thus, the nonlinear term in the corresponding GP equation, can be 
manipulated by varying the external magnetic field near the Fcshbach resonance [7] . This 
opens additional possibilities to control the quantum macroscopic dynamics of BECs. 

The NLS equation is a fundamental model for many physical media. A well-known ap- 
plication of the ID NLS equation is the description of the pulse propagation in nonlinear 
optical fibers [8]. Similar to optics, where bright and dark solitons are supported by the 
focusing and defocusing nonlinearity, respectively, in BECs the s-wave scattering inter- 
action between atoms is a similar determining factor. Therefore, condensates constrained 
to the one-dimensional shape can form dark or bright solitons. Dark solitons were found 
in condensates with repulsive interactions [9,10,11,12], while condensates with attractive 
interaction were recently experimentally shown to form stable bright solitons [13]. The 
appearance of solitons is one of the most interesting manifestations of nonlinear dynamics, 
and in the case of BECs it is of paramount interest, as in this case the solitons represent 
self- localized "waves of matter" . 

The similarity of the NLS and GP equations suggests that many nonlinear phenomena in 
BECs may have their counterparts in nonlinear optics, in particular, in fibers. However, 
due to unique manageability of the properties of BECs, new setups can be studied, which 
were not realized in nonlinear optics. One of such setups is related to the above-mentioned 
possibility of the effective control of the scattering length in BECs, i.e., the sign of the 
nonlinearity in the governing equations. 

Anticipating such experiments, in the present paper we find and study, analytically and 
numerically, a new family of solitons in two weakly coupled effectively ID condensates, 
trapped in a double-well magnetic potential, when the s-scattering lengths of the two 
condensates have different magnitudes, with particular emphasis on the case when the 
scattering lengths have opposite signs. 

The spatially separated condensates can be created by focusing a far-off-resonant intense 
laser beam, which generates a repulsive optical dipole force, into the center of a magnetic 
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trap. We assume that the condensates have the cigar-hkc shapes, with the transversal 
dimensions being strongly constrained by the trap, see Fig. 1. The chemical-potential 
difference between the traps, /iq, can be managed, moving the position of the barrier- 
generating laser beam by means of electro-optic or acousto-optic modulators (see also 
Ref. [14]). Such a setup was reahzed in experiments (see, for instance, Ref. [4]). We 
neglect a variation of the trap potential along the longitudinal (say, x) direction assuming 
it to be weak, so that for localized solutions well inside the trap the external potential 
can be considered as flat. Finally, as we are interested in the soliton solutions, we take 
into account the kinetic energy contribution. 

The corresponding coupled-mode equations, after obvious scaling transformations, may 
be cast in the following dimensionless form 

idtu + d^u + \u\'^u + V = 0, (1) 




where //q accounts for a difference of the chemical potentials between the two condensates, 

described by the order parameters u and v. It is assumed that the nonlinear interaction is 
always attractive in the -u-condensate, while in the t'-condcnsate the strength and sign of 
the nonlinearity are controlled by the coefficient a (e.g., repulsive nonlinearity if a > 0). 
The derivation of the couple-mode equations similar to Eqs. (l)-(2) from the correspond- 
ing GP equation was discussed in many works (consult, for instance, Refs. [15,16,17]), 
the outline is placed in appendix A. Thus, the system (l)-(2) realizes an interesting in- 
terplay between the dispersion (i.e., the kinetic energy contribution), self- focusing and 
self-defocusing nonlinearities, linear coupling and the potential shift, which deserves the 
study. 

Dynamics of two BECs in a magnetic trap was studied before. For instance, existence 
of a macroscopic quantum-phase difference, experimentally demonstrated in Ref. [4] , was 
used to study the coherent atomic tunnelling between two weakly coupled BECs confined 
in a double- well potential [16,17]. In Ref. [15] it was shown that the coherent oscillations 
due to tunnelling are suppressed when the number of atoms exceeds a critical value. The 
effect of the trap oscillations on the atomic tunnelling between two BECs was analyzed 
in Ref. [18]. In Ref. [19] an analogy with the nonlinear-optical directional fiber couplers 
was used to account for the kinetic terms in governing coupled NLS/GP equations (for 
a review of nonlinear dynamics in optical fiber couplers see section 6 of Ref. [20]; this 
analogy will play an important in the present work too). Further use of the analogy with 
guided-wave optics led to the development of a nonlinear collective-mode theory for BECs 
trapped in an external potential [21]. In Ref. [22], nonlinear modes in the form of chains of 
bright and dark solitons, which have no linear counterparts, were shown to be stationary 
solutions to the NLS equation with a multi-well external potential. A related model, based 
on a multi-component (i.e., vector) NLS equation, that describes repulsively interacting 
two-component BECs, possesses a coupled dark- bright soliton solution [23] . 

For a < (i.e., when the nonlinearities are focusing in both subsystems) equations sim- 
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ilar to the system (l)-(2) have been studied in connection with the nonhnear dynam- 
ics in the dual-core optical fibers (see the above-mentioned review [20] and original pa- 
pers [24,25,26,27,28,29]). However, the system (l)-(2) with a > 0, i.e., with opposite signs 
of the nonhnearity in the two condensates, was not considered before. Besides its rele- 
vance to the description of the coupled BECs, as argued above, the study of the model 
with positive a is of general interest: as we demonstrate in this work, it gives rise to 
interesting soliton solutions which exhibit some imusual properties (see figures 4 and 6 
below, for instance). Whereas letting the parameter a to have arbitrary negative values 
we also uncover novel features of the sohtons and bifurcations studied before for the case 
of a = —1. 

The paper is organized as follows. In the next section, we aim to identify a domain of 
the soliton existence (mainly for the case of a > 0) by deriving equations and inequalities 
which must be satisfied by the soliton solutions. In this section, we also find particular 
(sech-type) soliton solutions for a < 0. Generic numerically found soliton solutions are 
presented in section 3. In particular, a noteworthy result reported in this section is the 
occurrence of an unusual bifurcation (which was found for a = 1): a soliton may have 
its amplitude vanishing and width simultaneously diverging, while the number of atoms 
in this configuration remains finite. In section 4, we derive a criterion of the Vakhitov- 
Kolokolov (VK) type for the sohton stability, and study the fate of unstable solitons 
subject to small perturbations by means of direct numerical simulations. As the result, 
it is found that, depending on the values of the system parameters and on the form of 
the initial small perturbation, the unstable soliton either gives rise to a breather (with 
persistent intrinsic vibrations whose amplitude is dependent on the slope of the number 
of trapped atoms vs. the chemical potential) or completely decays into radiation. The 
concluding section summarizes results obtained in the work. 



2 The region of soliton existence and analytical solutions 

To search for the stationary sohtary-pulse solutions to Eqs. (1) and (2) we set 



where jJL is the normalized (dimensionless) chemical potential, and U{x) and V{x) are 
real functions vanishing as x — ±cx3. Thus we arrive at a system of two real ordinary 
differential equations for the pulse profiles: 



u{t, x) = e-'^^U{x), v{t, x) = e-'^^Vix), 



(3) 



+ {li + U^)U + V = Q, 



(4) 



dx^ 



+ {li-lio- aV^)V + [/ = 0. 



(5) 



dx^ 
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Equations (4)-(5) were numerically solved, to look for solitons in a wide domain in the 
parameter space (a, ^q, fi). However, before proceeding to the numerical solution, some 
results can be obtained in an analytical form, which will provide for some insight into the 
existence of solitons. The analytical results will make it possible to narrow a domain in 
the parameter space where it makes sense to look for solitons numerically. Moreover, the 
analytical results will be used to check the numerical solutions. 

First of all, considered as a dynamical system, equations (4)-(5) possess a Hamiltonian, 



Evidently, solutions vanishing as oo correspond to H = 0. We specify the class of 

the soliton solutions we look for as even solutions, U{—x) = U{x), V{—x) = V{x), with a 
single maximum at a; = (we aim to consider single- humped, i.e., fundamental solitons). 
For the definiteness' sake, we set U{x = 0) = Uq > 0. As the first derivatives of the fields 
vanish at the central point, we can derive the following quartic equation for the soliton 
amplitudes Uq and Vq = V{0) from Eqs. (4) and (5): 

[f^ + f ) U'o + (/^ - /^o - ^) V,' + 2UoVo = 0. (7) 



Further, multiphcation of Eq. (4) by dU/dx, and (5) by dV/dx, and the integration from 
X = to X = oo leads to the following identities: 



In principle, both in-phase (Vq > 0) as out-of-phase (Vq < 0) fundamental solitons are 
possible (recall we have set Uq > 0). Evidently, the right-hand sides of Eqs. (8) are negative 
for the in-phase, and positive for the out-of-phase sohtons. Thus, we conclude that the 
fundamental solitons must obey two inequalities: 



C/o'<-2//, if Vo>0, (9) 
aVo^<2{ii-iio), if Vb<0. (10) 

For a > (the opposite signs of the nonlinearities in the two subsystems), from these 
inequalities it follows that, for the in-phase fundamental solitons, the chemical potential 
of the [/-condensate is negative (// < 0), while, for the out-of-phase solitons, the chemical 
potential of the ^-condensate is positive (// — //q > 0). 

The domain where the soliton solutions are permitted is determined by the dispersion 
law. Assuming the exponential decay U ~ e~''^, V ~ e"''^ for \x\ —>■ oo, and linearizing 
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equations (4)-(5), we find two branches of the dispersion relation for the sohtons: 



(11) 



It is seen that the condition A;^ > imphes ji < jiQ, which excludes the out-of- phase sohtons 
for this branch for positive a. The appealing conclusion that the /ci-branch corresponds 
to the in-phase soliton solutions is confirmed by our numerical results (for both positive 
and negative a), see the next section. This fact is also used in the proof of the Vakhitov- 
Kolokolov stability criterion for the in-phase solitons in section 4. 

Finally, analysis of solutions of Eq. (7) by means of the inequalities (9) and (10) shows 
that, for fixed a > 0, /UQ) A^j and f/o, there is one negative [satisfying the condition (10)] 
and, at most, two positive [satisfying the condition (9)] real solutions for Vq. Thus, there 
may be at most two branches of the in-phase solitons (however, we were able to find, in a 
numerical form, only one in-phase soliton solution corresponding to a given Uq for fixed 
a, //o, and see the next section). 

For a < 0, when the system (4)-(5) describes two coupled condensates with attractive 
interactions, it also provides for a straightforward generalization of the model describing 
the (mismatched) nonlinear dual-core optical fiber. In this context, the soliton solutions 
were studied in detail for the case a — —1, i.e., equal Kerr coefficients in both cores 
[24,25,26,27,28,29]. In this case, the parameter jiQ accounts for the phase- velocity difference 
between the cores [27,28]. 

For negative a, Eqs. (4)- (5) have special exact soliton solutions, both in- and out-of- 
phase ones, which generalize, respectively, the well-known symmetric and antisymmetric 
solitons in the standard model of the dual-core optical fiber [24]. The exact solutions for 
the in-phase solitons are 



which exist in the special case, when and a are not independent parameters, but are 
related as follows: 




(12) 




(13) 



The out-of-phase sohtons are 




(14) 
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and they exist in the case 



(2) = 



1 



(15) 



'max — 



a 



[recall fio (a) is defined in Eq. (13)]. Note that the in-phase sech-type solitons correspond 
to the first branch of the dispersion law (11), Al/2 — k\, while the out-of-phase solitons 
correspond to the second branch, A\/2 — Unfortunately, the solutions (12) and (14) 
cannot be continued to positive values of a. 

Some additional exact results can be obtained concerning families of soliton solutions to 
Eqs. (4)-(5) with a < 0. For instance, by using the standard bifurcation analysis it is 
easy to show that there is another family of the in-phase solitons, bifurcating from the 
solutions (12) at the point 



if Hq and a are related as in equation (13) (this bifurcation is considered in detail nu- 
merically in section 3.2 and analytically in section 4). Setting a = —1, one recovers the 
bifurcation which was found, in terms of the dual-core fiber, in Ref. [24] (in this case 
/Xo = 0). Note that due to the relation /Xbif < A'-mL (which can be easily checked to hold) 
this bifurcation is present for any a < 0. Actually, the bifurcating in-phase solitons gener- 
alize the so-called A-type solitons (defined in Ref. [24] for a — —1, see also Refs. [27,28]) 
for arbitrary (negative) values of a. 

Another family of sohton solutions, the so-called 5-type sohtons, bifurcates from the out- 
of-phase solitons in a non-perturbative way [24] . However, the B-type solitons were shown 
to be unstable in the whole domain of their existence, while the out-of-phase solitons could 
be stable only in a very narrow interval [25]. For this reason, we discard the out-of-phase 
solitons from further consideration, so that all the solitons considered below are implied 
to be of the in-phase type. 



3 Soliton solutions (numerical analysis) 

We will discuss only the normalized quantities corresponding to the system (l)-(2). The 
numbers of the trapped atoms are also normalized accordingly, so that our discussion 
below is general instead of being tied to a particular shape and size of the external trap. 
The actual numbers of trapped atoms depend heavily on the particulars of the trap 
according to formula (A. 8) given in appendix A. The transformation coefficient between 
the physical and "normalized" numbers of atoms is estimated there to be on the order 
of 10^ for the current experimental setups. Thus the soliton solutions we discuss below 



4 - a 



(16) 
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correspond to the interval from thousands to a hundred thousands of the trapped atoms. 
Those are the characteristic numbers of atoms in the current experiments. On the other 
hand, it is not surprising to come up with precisely this interval in our model, since it is 
derived under the assumption of the weak nonlinearity, in which case the upper bound 
on the number of trapped atoms is determined by the external trap. To avoid confusion, 
we note here that below by the "number of trapped atoms" we mean only the number of 
particles in the model (l)-(2). 

To solve Eqs. (4)-(5) numerically the following iterative scheme was adopted 



^ ^ + ([/("-I)) + V^(-) = (17) 



dx^ 
dV^") 
dx^ 



+ a + = (18) 



i.e. we solve the eigenvalue problem for /x*^"^ at the n-th iterative step. Selecting the lowest 
eigenvalue at each iterative step leads to convergence to a sohton solution ^ (the use of 
this scheme, and the selection of the eigenvalue, are prompted by the well-known facts 
from the one-dimensional eigenvalue theory). We used the Fourier spectral (collocation) 
method with up to 256 grid points (see Refs. [30,31,32] for introduction into the spectral 
methods). Geometric convergence of the numerical solution to a soliton was noted for 
quite arbitrary initial profiles. As our method required only one of the two amplitudes 
to be specified (by appropriate normalization of the eigenfunction at each step), while 
the other was the result of the computation, we used the exact relations (7) to check the 
correctness of the numerical solutions. We also tested our approach, using the explicit 
soliton solutions (12). 



3. 1 Solitons for a > 



In the model with the opposite signs of the nonlinearity in the two condensates, a > 0, 
we have found a family of soliton solutions existing for all values of a and /iq, and all the 
values of the chemical potential fi permissible by the dispersion law, i.e., which satisfy the 
condition 




It should be stressed that, in comparison with the previous works done in the context of 
nonlinear optics, these are essentially novel solitons, as all the previously studied cases 

^ In some cases for a < 0, when there are several soliton solutions corresponding to the same a, 
jjLQ, and 11, wc used the shooting method to obtain the branches of solitons, to which the iterative 
scheme (17)-(18) has poor convergence. 
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[24,25,26,27,28,29] assumed the same sign of the nonhnearity in both cores of the system 
(ahhough a possibihty of the existence of bright gap sohtons was shown in the case when 
the signs were opposite in front of the second derivatives [29]; in optics, this case is quite 
possible, corresponding to opposite signs of the group-velocity dispersion in a dual-core 
fiber with asymmetric cores, while in the case of BEC this case makes no sense). 

Case I. a — 

We start with the special case a = 0, when the second condensate is characterized by the 
zero scattering length (no nonlinearity) . In figure 2 we plot the total number of atoms, 
and the numbers of atoms in each condensate as functions of the chemical potential /i for 
several values of the chemical-potential difference /iq. These dependencies, besides being 
important to quantity the BECs, determine the stability of the soliton solution: for fixed 
a and /iq (a > 0), the solitons are stable if the slope of the curve N = iV(/i) (A^ is the total 
number of atoms) is negative, dN/dn < 0. This is a stabihty condition of the Vakhitov- 
Kolokolov (VK) type [33] (proof for the case under consideration is given in section 4). 
At = yUinax, whcrc thc curvcs in Fig. 2 start from the zero number of atoms, the solitons 
bifurcate from the trivial solution U = V = 0. As /j, approaches //max, the width of the 
solitons tends to infinity, while their amplitudes decrease to zero. 

Decrease in the chemical potential difference //q causes the curve N — N{n) to sag in 
and develop a local maximum close to the upper limit value = yUmax, see Figs. 2(b)-(d) 
[in Fig. 2(d) we show only the sagging part of the curve, the shape of the whole curve 
being similar to that in Fig. 2(c)]. It is interesting to note that for n close to the value 
that corresponds to the maximum oi N — N{ii) almost all the atoms are trapped in the 
linear t'-condensate. The share of the number of atoms trapped in the condensate with 
no atomic interaction grows even further (for yU around the local maximum of = N{fj,)) 
with further decrease of /iq towards minus infinity. For instance, for a = and yUo = —10, 
the number of atoms in the linear condensate can be up to 99% (the respective function 
N — N{ijl) is similar to that in Fig. 2(d), but with the local maximum at A" 2500). On 
the contrary, for any given yUo, for large negative values of the chemical potential /i almost 
all the atoms are trapped in the li-condensate. 

The deformation of the solitons with the variation of the chemical potential for fixed a 
and yLto is illustrated by Fig. 3. The sohton solutions are plotted for three values of the 
chemical potential, which are marked by stars in Fig. 2(c), that correspond to the same 
total number of atoms. Note that, while at point 1, where ^ = —3.9894, the -u-component 
of the soliton is significantly higher, the f -component takes over as one moves to the right 
along the curve A" = A^(At) in Fig. 2(c): for yu, — —2.0236 and fi — —1.6576 (the points 2 
and 3, respectively) the T;-component of the soliton solution has a larger amplitude. 

II. Small positive a 

For small positive values of a (roughly, for a < 0.1) and positive or small negative /iq, 
the numbers of atoms in the two condensates vs. the chemical potential have the forms 
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similar to those of Figs. 2(a) through (c). For instance, for a = 0.1 and for no = —1, the 
curves are found to be similar to those in Fig. 2(c). Thus, for such values of a and /iq, 
the solitons bifurcate from the trivial solution at /x = /Xmax (19), as in the case of zero 
a. Moreover, for very small values of a and some negative fj^o, almost all the atoms can 
be trapped in the condensate with the (weak) repulsive inter-atomic interaction. As an 
example, a zoomed- in part of the functions = N{fi), Nu = Nu{fi) and = N^{fi) 
around the local maximum, together with the soliton solutions, are displayed in Fig. 4 for 
a = 0.001 and /iq = —5. At point 2 in this figure, about 96% of all the atoms are trapped 
in the f -condensate with the weak repulsive interaction [a zoomed-in part of Fig. 2(d) 
would be similar to left panel of Fig. 4 with a slightly lower local maximum] . 

However, the positive scattering length (self-repulsion) in the f-condensate brings new 
features too, as compared to the linear w-condensate. Most importantly, the grows of the 
local maximum of A" = A^(/x) with decrease of eventually saturates and changes for 
decrease, i.e., the maximum reaches its peak value at some negative /iq. Further decrease of 
/iq causes the local maximum to disappear. Instead, a divergence develops at the boundary: 
dN/dn — s> oo as /X — >■ /Xmax- For instance, for a = 0.1 and fiQ = —2 there is no local maxima 
at all, while the share of atoms in the f-condensate is large for /x close to /Umax due to the 
above mentioned divergence. Recall that the VK stability criterion demands the negative 
slope of dN/d/i, hence the solitons with a large share of atoms in the repulsive condensate 
are unstable for such values of a and /xq- This case is illustrated by Fig. 5, where we give the 
numbers of atoms vs. chemical potential and typical u- and v-shapes of the corresponding 
unstable soliton with larger number of atoms in the repulsive condensate. 

The following property of the soliton solutions is observed as /i — > fJ^max {ioi curves similar 
to those in the left panel of figure 5): the soliton's amplitude remains bounded, while its 
width is not. Therefore, our conjecture is that the solitons develop a "pedestal" (long 
shelf) of an increasing width, as the chemical potential approaches its limit value /Xmax 
(this limit is very difficult to study numerically, precisely for the same reason). In any 
case, all such solitons are unstable, since the slope of the corresponding curve A^ = A"(/x) 
is positive, see Fig. 5. 

We have checked that, for various small positive a and sufficiently large negative values of 
/xo, the dependence of number of atoms on the chemical potential is quite similar to what 
is shown in Fig. 5. Vice versa, for arbitrary negative /iq there is a threshold value of a such 
that, for a larger than this value, the derivative dN/d/j grows to infinity as — * /Xmax- 
For example, for /xq = — 1 it was found that, for a = 0.25, the dependence of the numbers 
of atoms on the chemical potential has essentially the same form as in Fig. 5. 

The results on the solitons for small (positive) values of a and various /Xq can be summa- 
rized as follows: the soliton solutions with fi close to /Xmax and positive or small negative 
fiQ are similar to those displayed in Fig. 3, and for large negative /xq are similar to the 
soliton displayed in the right panel of Fig. 5. On the other hand, for /x -C /Xmax the solitons 
are similar to the solution with /x = —3.9894 displayed in Fig. 3 (tagged by 1). 
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III. Large positive a 



For a ~ 1 and negative values of /iq, the curves N = iV(/x), iV„ = A^u(/x), and A^^ = Ny{ii) 
are still similar to those in Fig. 5, while for positive /Iq they are similar to the curves 
in Fig. 2(a). Increasing a further (i.e., making the repulsion between atoms in the v- 
condensate still stronger), we have found that the numbers of atoms in the condensates 
vs. the chemical potential take the form of what is displayed in the left panel of Fig. 5 
also for positive iiq. 

The value separating the two different types of behavior of iV = N{fi) as // ^ /imax- i-c, 
given by Fig. 2(a) and the left panel of Fig. 5, was found to be a = 1. At a = 1 and 
/^o = 0, the soliton has its amplitude gradually vanishing, as /i —>■ but, nevertheless, 
the number of atoms approaches a non-zero value, which is evident from the upper panel of 
Fig. 6. This seemingly strange bifurcation implies that in the limit /i — > //max the soliton 's 
width is coupled to its amplitude, thus the integral that gives the number of atoms remains 
constant (see the bottom panel in Fig. 6). Such a bifurcation may be naturally called a 
discontinuous one, as the soliton bifurcates from the zero solution without a visible jump 
in the shape, but with a jump in the number of atoms. 

It may be interesting to check if the numerically found solitons for a > are approximated 
by a sech-based ansatz. Let us assume, for example, that the soliton solutions can be 
approximated by the sech" -functions, i.e., the solitons have the shape of / = Asech.^ [x / d) , 
where A (amplitude) and d (width), in general, different for the two components, are the 
dynamic parameters (determined through the variational analysis) and the power a is 
fixed. This is an improvement of the usual sech-type ansatz. Given a numerical soliton 
solution we need to determine a such that the corresponding sech°-function approximates 
the soliton in some sense. This can be done in many different ways. Here we adopt as the 
proximity measure the following functional 

oo oo 

Pf = 5^^^ • (20) 

J dx f J dxx'^p 

— oo — oo 

Inverting the function P — Pgech" = P{ot) and using the numerically computed Ps (for a 
numerically found soliton solution S — S{x)) we get the corresponding a: a — a{Ps)- [The 
proximity measure was chosen in the ad hoc way. However, there arc two conditions to be 
satisfied: (i) the functional must be independent of both A and d; {ii) the function P = 
P{a) must be monotonous (it is decreasing for Pf).] If the solitons could be approximated 
by the above ansatz, the expected result would be a sharply peaked distribution of the 
power parameter values a computed for the numerically found solitons. Then one could 
assume the value of a corresponding to the peak as the approximation value. 

We have computed the power parameter a in the above approximation for the numerically 
found soliton solutions for various a, /iq and /i. For instance, for the solitons shown in 
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Fig. 3 the corresponding values of a are as follows (the subscripts refer to the u- and 
v-components of the soliton) 

= 0.825, = 1.61, 
Q;(f) = 0.218, q;(2) = 1.01, 
a^f) = 0.711, q;(^) = 0.922. 

Evidently the soliton solutions cannot be approximated by the sech"-type functions with 
the same common power a due to a broad distribution of the power parameters corre- 
sponding to the solitons. Here we note that, in contrast, the variational approach based 
on the sech-type ansatz with the amplitudes and widths as the dynamic variables was 
successfully applied before to the solitons for a < (see Refs. [27,28,29] where the case of 
a — —l was considered). 

3.2 Solitons for a < Q 

In the case of a < 0, when the interactions are attractive in both condensates, it is enough 
to consider the interval — 1 < a < 0, as Eqs. (4)-(5) are invariant against the following 
substitution: 

U ^ U = \J —aV^ V ^ V = \/—aU, Uo ^ P'O = ~l^o, a —>■ d = -, (21) 

a 

and II ^ jl — II — Ho- 

The solitons for a < were studied before in the context of the nonlinear dual-core optical 
fibers [24,25,26,27,28,29], chiefly in the particular case a — —1. Though our new results, 
presented below, pertain to a ^ —1, we also consider in detail the previously studied case 
of a = — 1 and obtain some results, which were not known before. 

We will frequently refer to the special case when a and //q are related as follows /Xq = 
-\/^ — = Ho (a). Note that only in this case Eqs. (4)-(5) admit (real) sech-type 

soliton solutions. 

In section 2 it was shown that there are two branches of the in-phase solitons, one given 
by Eq. (12) for arbitrary (negative) a and //q = A*o(ct)) and the other one bifurcating 
from these sech-type solitons at iJ, — //bif, see Eq. (16). We have found numerically that 
these branches coexist for arbitrary values of a and iiq, although they do not always 
intersect, hence do not always undergo a collision bifurcation. Thus, the sech-type soliton 
(12) is a special case of a broader family of soliton solutions. Moreover, we have verified 
numerically that this broad family of (in-phase) solitons corresponds to the /ci-branch of 
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the dispersion relation, see Eq. (11), i.e. the sohtons exist for /i satisfying the condition 
(19). 

To present the results in more detail, we first consider small (negative) values of a. We note 
that the bifurcation involving the two branches of solitons, with one branch corresponding 
to the sech-type solitons (12), belongs to the tangential type, i.e., it is a one-sided cusp 
bifurcation (consult, for instance, Ref. [35]). This is clearly seen from Fig. 7, where we 
take the values a = —0.1 and /xq = Aio(ct) (this statement is also proven analytically in 
section 4). For small deviations /xq — IJ'o{o,), and small a (roughly for —0.5 < a < 0), the 
bifurcation is similar. For larger deviations /Iq — iJ>o{(i), the two branches develop a trend 
to collide, which gives rise to a picture resembling a collision bifurcation (see Fig. 8, where 
this type of behavior is shown, though for a larger negative value of a). 

Finally, for small negative a and for small but finite deviations /Xq — A*o(a) we have checked 
that the solitons belonging to the branch which corresponds to the sech-type sohtons at 
A^o = /^o(fl) can be uniformly approximated by the sech ansatz for all values of the chemical 
potential /x. This is, however, not so for larger negative values of a (see the discussion below 
for a — —1). 

We have found that the collision bifurcation is replaced by the coexisting branches when 
the deviation /xq — /xo(a) has large enough modulus. For instance, when a = —0.8, the two 
branches are clearly separated for fiQ = —0.2, see Fig. 8, which is ~ 10% away from the 
corresponding value ^^{a = —0.8) ~ —0.22 (the case of the small deviation /iq — /xo(«) is 
analyzed in detail for a = —1 below). 

In Fig. 9 we display the bifurcation diagrams for a = —1. Here it is worth noting that, 
due to the symmetry against the substitution (21) the solution components U = U{x) 
and V = V{x) are interchangeable in this case. For instance, the sech-type solitons have 
identical components, U = V (the symmetric solitons), while the solitons belonging to 
the other branch (the asymmetric ones) admit two possibihties: U > V, or V > U. Only 
the diagrams with U > V are displayed in Figs. 9(b) and (d). 

First of all, we notice that Figs. 9(a) and 9(b) reproduce the known bifurcation diagram 
for Ho = [24] (note that in this case iiq = //.o(a = —1)). Although it is demonstrated in 
section 4 that the asymmetric solitons bifurcate tangentially from the symmetric branch, 
we were not able to show this feature in the picture, because a region where it takes 
place for a = —1 and //q = is extremely narrow (for the same reason, this fact was left 
unnoticed in Ref. [24]). 

As soon as the chemical-potential difference /Iq deviates from zero (in terms of the dual- 
core- fiber model, this corresponds to a phase mismatch between the cores), the collision 

bifurcation undergoes a significant transformation, see figure 9(c). It is interesting to note 
that, as soon as /iq goes away from zero (we take /iq = —0.01), the former sech-type branch 
divides itself into two different branches of soliton solutions, while the former asymmetric 
branch splits into two close curves, which pertain to different branches as well. 
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Consider the top part of the branch with the arrow in Fig. 9(c). In this region, the power 
parameters and «„, numerically computed using Eq. (20), differ from I by less than 
2% . Although this curve corresponds to the solitons which are well approximated by the 
sech ansatz and the corresponding soliton solutions have nearly identical components, 
[/ ~ y, it nevertheless can be identified with the deformed asymmetric branch of solitons. 
This statement is supported by comparison of Fig. 9(d) with Fig. 9(b): in both cases, the 
bifurcation diagram has a pitchfork- type form, but in the case of Fig. 9(d) the stem and 
the middle prong of the pitchfork, both lying very close to the sech-type branch for hq = 0, 
belong to different branches of solutions. We were not able to trace the bifurcation point 
A*bif by continuing Fig. 9(c) in the direction of negative /i, i.e., for /i — > — oo, and checking 
stability of the solitons (the bifurcation point is the threshold of the soliton instability, 
see section 4). Thus, when a = — 1, the collision bifurcation is replaced by two coexisting 
branches of solitons already for hq = —0.01. 

Finally, we observe that the branch of the sohton solutions denoted by the thin curves in 

Figs. 7 through 9, which corresponds to the asymmetric solitons for a = —1 and /xq = 0, 
undergoes a turning-point bifurcation. The turning point is clearly seen in these figures. 

As it was mentioned above, we would not consider bifurcations of the out-of-phase (e.g., 
antisymmetric) solutions, as they are always unstable (see also Ref. [25]). For the same 
reason, we do not consider of higher-order multi-humped solitons, which are also found 
in the numerical solution, but turn out to be unstable in all the cases. 



4 Soliton stability and evolution under perturbations 

After having found the soliton solutions, the next necessary step is to analyze the soliton 
stability under action of perturbations, as well as the fate of unstable solitons. Accordingly, 
we split this section into two parts. In the first part, we prove the above-mentioned VK 
stabihty criterion, dN/d/j, < 0, for the in-phase soliton solutions of Eqs. (4)-(5) with 
a > and show where the VK criterion applies to the solitons for a < 0. For instance, we 
demonstrate that the bifurcation point (16), discussed in the previous section, is (quite 
naturally) an instability threshold for the sech-type solitons (12). Our proof of the stability 
criterion for the soliton solutions of the system (l)-(2) generahzes a similar proof for the 
solitons of a single (scalar) equation, which can be found in Refs. [33,34]. In the second 
part of this section, we numerically study the evolution of unstable solitons subject to 
initial perturbations. 

Here we point that for a — —1 (with hq — 0) the sech-type soliton solutions (12) discussed 
in the present paper, and the in-phase solitons in general, reduce to the solitons which 
were studied in Ref. [24] , in the context of the dual-core optical fibers. The stability prop- 
erties of the in-phase solitons for a = —1 were investigated in Ref. [25], and evolution of 
unstable solitons under perturbations was simulated in Ref. [26]. Below, we (in particular) 
generahze the stability results of Ref. [25] for arbitrary a < and //q, which is relevant to 
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the model (l)-(2) for trapped BECs. 



4-1 Stability analysis 

I. Stability of the in-phase solitons for a> 

The sohton stabihty with respect to small perturbations is determined by the system 
(l)-(2) linearized about the soliton solution. Thus we need to consider a perturbed soliton 
solution: 

u = e-^^* [U{x) + ui{t, x)], e-^^* [V{x) + Vi{t, x)] , 

where U{x) and V{x) are the stationary soliton solution components, while Ui{t,x) = 
UR{t, x) + iUi{t, x) and vi{t, x) = Vi?(t, x) + iVi{t, x) represent the perturbation. Lineariz- 
ing equations (l)-(2) with respect to the perturbation we arrive at the system 







lJ") 


-l\ 




Ui 




-lS") 


1 


Ui 


Vr 




-1 


L^^ 


Vr 




V 


1 


-lS") j 





(22) 



where the linear operators are defined as follows: 

4"^ ^ - + + ^U\x)^ , lS^) = - + - /.o - ?>aV\x)^ . (23) 

The solitons are unstable if the matrix operator on the right-hand side of equation (22) 
has an eigenvalue A with the positive real part [in fact, the eigenvalues are real, since 
Ao is non-negative, see Eqs. (24)-(25) and the discussion below]. The system (22) can be 
written in an equivalent second-order form, hence the corresponding eigenvalue problem 
can be reformulated for a fourth-order differential operator acting on a two-component 
vector (X„,X^), with and defined as follows (consult also appendix B) 

UR{t, x) = e-'^^Xu{x) + c.c, VR{t, x) = e-'^^X^{x) + c.c. 
Here, — i\ has the meaning of a frequency. The linear eigenvalue problem then takes 
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the form 



AoAi 



(24) 



where we have introduced two symmetric matrix operators: 



Ao 



^0 

-1 h 



{V) 





Ai 



LP -1 



^ -1 L\ 



(25) 



The proof of the stability criterion rehes on properties of the factorization operators Aq 
and Ai. We first aim to show that Aq is non-negative, i.e., the scalar product (^'|Ao|^'), 
defined as Xll"^^ dx^'^(a;)Ao^(x) (here | stands for the Hermitian conjugation), is non- 
negative for any real two-component vector ^{x) — {^\{x\'\\}2{x)Y . Indeed, the following 
inequalities follow from the definitions the operators Lq"^ and Lg^^ : 

/ Axi}\x)Lt\{x) > [ dx^^\i/j(x)\\ [ dxilj*(x)L^o'>ilj(x) > [ dx^^\i/j(x)\^ (26) 
J J Uix) J J Vix) 



To arrive at Eq. (26) it is enough to note that, due to Eqs. (4)-(5), these operators can 
be rewritten as 



■in) 



1 d 



U\x) 



d 1 



+ 



V{x) 



-{v) 



1 d 



U{x)dx ' 'dxU{x) U{x) 



'I ^0 



V\x) 



d 1 



+ 



U{x) 



V{x)dx ' 'dxV{x) V{x) 



Taking advantage of the inequalities (26) and of the fact that, for the in-phase soliton 
solutions, the product U{x)V{x) is positive, we get 

{*|A(,I*> = [dx {v.;(i)4"Vi(i) + V'2*(a:)4"'V'2(i) - riix)^^) - '/'i(i)'/'2*(i)} 



>/dx{ 



V(x) I , . . ,2 U(x) . , , ,.0 



I 



dx 



V{x) 



x) 



U{x) 



x) 



> 0. 



The operator Aq has, in fact, a zero eigenvalue, with the soliton solution proper ([/, V)^ 
being the corresponding eigenf unction: 



Ao 



U{x) 



\V{x) 



0. 



(27) 
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In its turn, the other factorization operator Ai also has a zero mode, 
\dV{x)/dx ) 



(28) 



which, however, has a node (zero) at a; = 0. Thus, according to the Sturm oscillation 
theorem, the operator Ai may have negative eigenvalues (it is found numerically that 
Ai has either one or two negative eigenvalues). In general, it may have two nodeless 
eigenfunctions corresponding to two negative eigenvalues: 



Ai 



/«( 





A, 




where f^^^g^^^ > and f^'^^g^'^^ < 0. However, it is shown in appendix C that Ai can- 
not have an eigenfunction of the second type for a > (corresponding to a negative 
eigenvalue) . 



Return now to the eigenvalue problem (24). We have 



AT,, 




(29) 



where Cq is a scalar constant. In deriving (29), we have taken into account that for 



{{u,v)m^o, 



(30) 



since Aq is a symmetric operator. Thus, the relation (29) is justified. The minimum eigen- 
frequency Qq can be also found from the following minimization problem: 



mm 



(^|Ai|^) 
(*|Ao'l*) ' 



(31) 



with ^ = {ipi,ip2)'^ subject to the constraint (30). To evaluate the sign of has 
to evaluate the sign of the numerator in the expression (31). Due to the constraint (30), 
the latter problem is equivalent to evaluation of the sign of the lowest eigenvalue of the 
following generalized eigenvalue problem. 



Ai 





= A 


f1 


+ /3 


f1 













(32) 
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where P is a. constant and ^ = {X, Y)'^ obeys Eq. (30). From the above consideration, we 
know that Ai has only one negative eigenvalue, and the cigcnfunction cannot be orthogonal 
to {U, V) for UV > 0, since f^^^g^^^ is positive too. Therefore, the expansion of {X, Y)^ 
over the eigenfunctions of Ai, 



Ai 



4") 



A,. 



4") 



has the following form (for simplicity of the presentation, the summation below also 
denotes the integration over the continuous spectrum) 




((0S"\0r^)i 




(33) 



where the zero eigenvalue A2 = does not enter the sum. In deriving Eq. (33), we have 
made use of Eq. (32) and the orthogonality of the eigenfunctions of Ai. The constraint 
(30) requires that the generahzed eigenvalue A from Eq. (32) satisfies 



2 



u 



An — A 



pF{X) = 0. 



(34) 



Consider the function F = F{^) defined in (34). It is discontinuous at the points { = A„ 
and is a monotonically growing function elsewhere. Note that the summation involves 
one negative eigenvalue Ai, the next (smallest) eigenvalue A3 being positive. Then, since 
the generalized eigenvalue A defined in Eq. (32) is a zero of the function F, we have 
sgnA = — sgn{F(0)}. Therefore, we must demand F{0) < for the soliton stability. Note 
that 



F{0) = {{U,V)\A^'\ \ ^ I). 



(35) 



On the other hand, the differentiation of Eqs. (4)-(5) with respect to n yields (note that 
the operator Ai is even with respect to the substitution x — > —x ) 



Ai 



^ dU/dA 
^dV/di,) 




(36) 
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or 



-(;:)=( 



dU/d^ 



Finally, substitution of the latter result into Eq. (35) leads to the VK criterion for the 
soliton stability: 

Au\ I dUldfi\ IdN 



It should be pointed out that, in the course of the derivation of the stability criterion (37) 
for the in-phase soliton solutions of Eqs. (4)- (5), we used the following two facts about 
the operator Ai: (i) there is only one (non-degenerate) negative eigenvalue, and {ii) the 
corresponding eigenfunction is not orthogonal to the soliton solution ([/, V) (otherwise the 
lowest generalized eigenvalue A coincides with the negative eigenvalue of Ai). With these 
two conditions satisfied, one can use the VK stability criterion for other in-phase soliton 
solutions of Eqs. (l)-(2). Below, this fact is used for the stability analysis of the sech-type 
solitons (12), and the solutions bifurcating from the sech-type solitons at /x = /Xbif- 



II. Stability analysis of the in-phase solitons for a < 



First of all, let us consider the special case when /xq is not an independent parameter, 
but is the function of a considered above, i.e. /iq = — l/\/^a = f^o{a), see Eq. (13). 
To derive the expression for the bifurcation point (16), we note that, in general, at such 
a point two branches of the soliton solutions {U^^\V'^^^) and {U^'^\V^'^^) coincide, hence 
due to Eq. (36) we have 



Ai 



( dU^^ydfi\ 
dV^''^/dii j 



(38) 



Therefore, the bifurcation point is characterized by the appearance of the second (node- 
less) zero mode of the operator Ai, in addition to the zero mode given by Eq. (28). Such 
a zero mode {Xi, X2)'^ can be easily found analytically for the sech-type solitons (12); in 
this case, it is a solution to the system 



d^Xi 
dx^ 



+ 



+ 



2„„„i,2 / 



V2 



+ ^2 = 0, 
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^2 + = 0. 
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Solving these equations, we obtain the expression (16) for the bifurcation point /luf, and 
the corresponding zero mode 



sech 



'2(1 



3^y—a 



1/2 



(39) 



Now, we aim to show that the sech- type sohtons (12) are stable in the interval 



A*bif < A* < A«; 



(1) 



(40) 



where n^l^ — —1/^/^ [see Eq. (13)], and unstable otherwise [we remind that here 
A*o = A*o(a)]- We need to determine the direction of the shift of the zero eigenvalue for 
7^ /^bif as we move along a given branch of the soliton solutions. To this end, we can 
use the perturbation theory for eigenvalues of linear operators. A shift of the chemical 
potential from the bifurcation point by a small amount e, fj, = fiuf + ^, results in a 
perturbation of the eigenfunction (39) and the corresponding eigenvalue A = 0+€li+O{e^) . 
Substitution of the perturbed eigenfunction. 




into the eigenvalue problem 
Ai ^ I = A 



^X, 



Xo 



Xo 



(41) 



and keeping only the first-order terms in e, we derive the following equation for li, 



Ai 



















[x,i 





1 + m{dU/d^) 

1 - QaV{dV/dii) ) 



\X2j 



(42) 



where the soliton solutions U and V are given by Eq. (12), and the operator Ai, together 
with V and their derivatives, are taken at = /^bif • The left multiplication of Eq. (42) 
by {Xi, X2) and integration over x leads to the following expression for the sign of li. 



sgnZi = -sgn y dx 1 1^1 + 6C/^j Xl + - ^aV^^ Xl 



(43) 



The expressions (42) and (43) are valid for the two branches of the soliton solutions which 
collide at the bifurcation point /^bif (16). 
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We now apply Eq. (43) to the sech-type solitons of Eq. (12). We have 
U — ^isechy, — — = — 2sech^|/(l — ytanhy), V — . , y = 



Aix 



V2' 



a 



thus 



sgn/i = — sgn 




6sech^7/[l 



7/tanhy])sech'^7/ > 0. 



Therefore, noticing that sgnA = sgn(e/i), we arrive at the following formula for the sign 
of the eigenvalue as we move along the sech-type branch of the soliton solutions: 



Thus, we have A > for /x > /Xbif- As the result, the nodeless eigenfunction from Eq. (41) 
with the property X1X2 < corresponds to a positive eigenvalue, and Ai has only one 
non-degenerate negative eigenvalue for //bif < jJ. < A*mL- 

To complete the proof of the stability of the sech-type solitons in the interval (40), we note 
that the eigenfunction corresponding to the (only) negative eigenvalue Ai of Ai cannot 
be orthogonal to the soliton solution {U,V) (12). Indeed, such an eigenfunction (^1,^2)^ 
has no nodes and satisfies the condition > at /i = /Xbif, as it is orthogonal to the 

zero mode (39). Thus, to become orthogonal to the soliton solution (f/, V), one of the 
components of this eigenfunction should first pass through zero at some value of fi which 
is impossible. 

For ^0 = Ato(fl)) the stability properties of the in-phase solitons bifurcating from the sech- 
type soliton solutions at = /Xbif arc affected by the turning-point bifurcation mentioned 
in section 3 (see, for instance, the top panel of figure 7). To find where the VK criterion 
applies to these solitons, we note that, at the collision bifurcation point //bif , these soliton 
solutions share the operator Ai with the sech-type solitons and at the turning point one 
of the (two) negative eigenvalues of Ai passes through zero and becomes positive, as we 
go from the upper part of the branch to the lower one [this is due to the fact that the 
function F(^) makes a jump from —00 to -|-oo at ^ = as we pass the turning point 
downwards, due to dN/dji — 2F(0)]. Therefore, as we move along this (i.e., bifurcating) 
branch from the bifurcation point /ibif, the zero eigenvalue corresponding to the nodeless 
zero mode (39) of Ai becomes negative, but at the turning point it passes through zero 
and becomes positive. Therefore, the VK criterion applies to the lower part of the branch, 
while the solitons of the upper part are always unstable (there is one negative generalized 
eigenvalue for them). 

To complete the consideration of the special case with /iq = /Lio(a), we note that, at the 
bifurcation point 11 = /luf, the numbers of atoms for the two branches of the in-phase 
solitons vs. chemical potential, say A^i = A^i(a*) and A^2 = A^2(a*), have equal slopes: 



sgnA|sech = sgn(/i - /Xbif)- 



(44) 
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dNi/djjL = dN2/dfj,, which is a straightforward corollary of Eq. (36). Indeed, the zero 
mode (39) of the operator Ai, also defined as 



is orthogonal to the soliton solution {U, V)'^ taken at the bifurcation point //bif , as it is 
the right-hand side of Eq. (36). Thus, 



Geometrically, it means that the in-phase solitons undergo a tangential (one-sided cusp- 
like) bifurcation. This fact was already illustrated numerically in the previous section. 

The stability properties of the soliton solutions for /iq ^ /xo(a) can be summarized as 
follows. According to their presentation in Figs. 7, 8 and 9(c)-(d), we will refer to the two 
branches of the soliton solutions as the thin and the thick branch, respectively, where the 
latter one corresponds to the sech-type solitons (12) for /iq — //o(a)- 

First of all, the VK criterion apphes to the sohtons belonging to the thick branch for all 
values of /i, thus they are stable when the number of atoms decreases with increase of 
the chemical potential. Second, to understand the stability of the solitons belonging to 
the thin branch, we note that, at the turning point, one of the negative eigenvalues of 
Ai passes through zero as one goes from the upper part of the branch to the lower part, 
similar as in the special case Ho — iJ,o{a). Thus, the solitons belonging to the upper branch 
are unstable, as the operator Ai has two negative eigenvalues. The VK criterion applies 
to the lower branch if the (only) negative eigenvalue of Ai there has the eigenfunction 
which is non-orthogonal to the soliton solution proper. It is enough to verify this condition 
numerically just at one point. The outcome is that the VK criterion is indeed applicable 
to the lower part of the thin branch. 

The analytical results on the soliton stability were checked by direct numerical solution 
of the linear eigenvalue problem (24), and by counting the negative eigenvalues of the 
operator Ai (we used the Fourier spectral discretization method with up to 256 grid 
points in the LAPACK routines of Matlab). As the outcome, it has been verified that the 
analytical results completely agree with numerical ones in all the cases. 

4-2 Evolution of perturbed unstable solitons 

Here we report results of direct numerical simulations of Eqs. (l)-(2) with perturbed 
soliton solutions as the initial conditions. We have used the Fourier spectral discretiza- 
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tion method in the x coordinate combined with the leap-frog time-stepping scheme. The 
stabihty of the scheme itself was guaranteed by selecting a small enough time step (we 
had At = 0.001), such that the stability domain of the leap-frog time-stepping method 
contains the eigenvalues of the spacial discretization operator multiplied by At. The ra- 
diation was absorbed by introducing a smooth distributed damping (smooth damping is 
necessary for stability of a numerical scheme based on the spectral methods), which had 
a negligible effect on the localized solution. 

Previously, the results of numerical simulations for the case of negative a in the context 
of the nonlinear dual- fiber optical model were reported in Ref. [26], where a — —1 and 
/xo = were used. Here, our main interest is the evolution of unstable solitons in the most 
interesting case, when a > 0, corresponding to nonlinearities of the opposite signs in the 
two condensates. 

First of all, we note that the in- phase solitons may have only one unstable mode, i.e., the 
linear eigenvalue problem (24) can have only one negative eigenvalue fl^ (or equivalently, 
one imaginary eigenfrequency fl). For a — —1 and /iq — 0, similar result was reported in 
Ref. [25]. 

In the context of the stability analysis presented above, it is easy to understand why 
there is just one unstable mode. Indeed, where the VK criterion applies, and the solitons 
become unstable due to the positive slope, dN/d/j > 0, this means that one generalized 
eigenvalue of the problem (32) has passed zero and become negative with the change of 
sign of -F(O), see Eq. (37). For the soliton solutions corresponding to the top part of the 
thin branches in Figs. 7-9 the VK criterion does not apply, but note, however, that the 
slope is negative. Thus, the only possible unstable mode is due to the appearance of the 
second negative eigenvalue in the spectrum of Ai and, hence, a single negative generalized 
eigenvalue of the linear problem (32). The fate of the corresponding (second) negative 
eigenvalue of Ai in this case depends on whether there is a collision bifurcation or not. In 
the former case, the negative eigenvalue must go to zero as we approach the bifurcation 
point /ibif- sec Fig. 7, since at /i = /^bif this eigenvalue is zero (the stability criterion 
applies to the branch of the sech-type solitons on one side from the bifurcation point). 
In the latter case, when there is no collision bifurcation, the negative eigenvalue tends to 
— cxD as we move along the top part of the thin branch away from the turning point, see 
for instance, the top part of figure 8. In both cases, this negative eigenvalue also goes to 
zero as we approach the turning point (i.e., as we move in the opposite direction). 

The fact that there is just one unstable mode, i.e., one unstable direction for the growth of 
a small perturbation, simplifies the task of understanding the fate of the unstable soliton 
solutions subject to perturbations: it is sufficient to solve the initial-value problem for 
Eqs. (l)-(2), taking as the initial condition the soliton plus a perturbation proportional 
to its unstable mode. 

Here we note that though the perturbation based on the unstable mode is complex (consult 
appendix B for more details), it is sufficient to use just its real part. This is due to the fact 
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that the sohton solutions are given by the real functions and, hence, only the real part of 
the perturbation, e.g., proportional to (Xu-X^)'^ from equation (24), affects the number 
of atoms in the first order approximation. Indeed, for a perturbation u — U = e(X„ + 
and V — V — e{Xy + iY^), where U and V are the soliton components, the corresponding 
variation of the number of atoms reads 

oo 

AN^2e J dx (UX^ + VX^) + ^(e^). (45) 

—oo 

Below by the small perturbation (or the unstable mode) we mean the real part. 

We now focus on the case of a > 0. We found that the evolution of an unstable soliton is 
determined by the shape of the stability curve N = N{ii) and the sign of the small pertur- 
bation in the form of the unstable mode. There are two distinct cases, which correspond 
to the shapes shown in Figs. 2(c) and the left panel of Fig. 5, respectively. 

We first consider the evolution of the unstable solitons in the former case. Figure 10 shows 
the (unstable) sohton and its unstable mode for this case. We chose a = (no nonlinearity 
in the t'-subsystem) and /^o = but the evolution is similar for other values of these 
parameters, provided that they give rise to a similar shape of the curve N — N{ii). For 
example, it is so for a and /iq of Fig. 4 (and in any other case of a > 0, when the v- 
condensate has atoms with repulsion, but the shape of N{fj,) is similar to the one given 
by Fig. 2(c)). 

Figure 11 illustrates the evolution of the unstable sohton solution, where the top picture 
indicates the direction of the evolution if one adds (the case denoted by 1) or subtracts 
(2) the unstable mode, as it was specified above. For instance, in scenario 1 there were, 
initially, more atoms trapped in the w-condensate, but this distribution of atoms in the 
condensates is unstable. For such initial conditions, the corresponding attracting config- 
uration has more atoms in the I'-condensate, as the arrow 1 indicates in the top part of 
Fig. 11. 

The picture in the bottom part of Fig. 11 shows the evolution of the two numbers of atoms, 
initiated by the instability of the original soliton. As it follows from this picture, in both 
cases the instability gives rise to a breather-like state with persistent internal vibrations. 
However, the amplitude of the vibrations is dependent on the slope of the stability curve: 
for the steeper slope the amplitude is smaller. 

Now let us consider what happens to the unstable soliton when the stability curve N — 
N{ii) has the shape of that shown in the left panel of Fig. 5. In this case, a scenario of 
the type 2 in terms of Fig. 11 (with formation of a breather) is also observed. However, a 
scenario of the type 1 does not take place in this case. It is replaced by complete decay of 
the soliton into radiation, as is shown in Fig. 12, which pertains to the case a = 0.25 and 
Ho — —1, the curve iV(/x) being in this case similar to the one in the left panel of Fig. 5. 
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5 Conclusion 



We have found novel soliton states for the system of two linearly coupled NLS equations, 
which describe two Bose-Einstein condensates trapped in an asymmetric double-well po- 
tential, in particular, when the atomic interactions have opposite signs in the two conden- 
sates. This system realizes an interplay between dispersion (corresponding to the kinetic 
energy in the condensates), attractive and repulsive nonlinearities, linear coupling, and 
the chemical-potential difference between the two traps. 

We have found stable soliton solutions with almost all atoms (~ 96%) being trapped in the 
condensate with weak repulsive interaction, which is a novel self-trapping phenomenon in 
a system of two weakly coupled condensates having the scattering lengths of the opposite 
signs. An unusual bifurcation was found too, when a soliton merges with the zero back- 
ground, having its amplitude vanishing and width diverging, while the number of atoms 
trapped in the soliton remains finite. 

The stability of solitons was studied in detail, and evolution of the unstable ones was 
investigated by means of direct numerical simulations. The outcome is that the unstable 

solitons either give rise to breathers or completely decay into incoherent waves (radiation). 
Our results for the two-component solitons are also relevant to the study of coherent 
atomic tunnelling (see, for instance, Refs.[16,17]), which has attracted a lot of attention 
in the context of the Bose-Einstein condensation. 

New results concerning the two-component solitons and their stability were also obtained 
for the asymmetric model in which the self-interaction is attractive in both components. 
This model has direct applications to the dual-core nonlinear optical fibers. Although it 
was studied in many works, new results for this model have been obtained here. 

The investigation of the solitons in linearly coupled quasi-one- dimensional condensates in 
a double-well potential can be continued in several directions. These include the study 
of breathers, and collisions between moving solitons. In particular, boundary effects will 
have to be taken into account for moving solitons, i.e., the form of the longitudinal po- 
tential will enter the stage. Moreover, the antisymmetric (out-of-phase) and higher-order 
(multi-humped) solitons, although unstable, may play an important role in the quantum 
tunnelling phenomena, since the solution may oscillate about these solitons during a finite 
time. These aspects of the soliton dynamics will be addressed elsewhere. 
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A Derivation of the coupled- mode equations 



Let us outline the derivation of the coupled-mode equations (l)-(2). We will use the tilde 

for the physical variables in the GP equation to distinguish them from the corresponding 
dimensionless variables of the system (l)-(2). The crucial observation is that the modulus 
of the order parameter (single-atom wave function) is exponentially small in the barrier 
region, see figure 1 of section 1. This allows us to approximate the solution to the cor- 
responding GP equation with the double- well potential 14xt(p) strongly confining in the 
transverse directions p = {y, z) (we neglect the longitudinal variation of Kxt), 

as a sum of the factorized wave functions for the two wells: 

^{i, x,p)^ x)^u{p) + Mi, x)^v{p), (A.2) 



where $u(p) and ^v{p) would be the ground-states (for the transversal degrees of freedom) 
if the two wells were isolated. Here we point that go is different in the two wells: g'o^Vfl'o"^ 
has arbitrary value and sign. This quotient can be easily managed by technique based 
on the Feshbach resonance [7] with the magnetic field applied to only one of the wells. 
Inserting the ansatz (A.2) in the GP equation (A.l) and using the conditions 

J ^M'pc^ 0, J $^dV= / ^^dV= 1, (A.3) 

one arrives at the system of two equations for the wave functions describing the longitu- 
dinal evolution of the condensates in the wells: 



ih-^ = -^-Q^ + i^v + gvlH^yipv - Ki^u- (A.5) 

Here the parameters are defined as follows (by changing the sign of either or <l>^, if 
necessary, one can set > 0) 

^u,v = j + ^e-t^i^) dV, 9u,v = J 5o<^dV, 

K^-J (^|^(Vp-$„)(V,-$,) + $„Kxt$.) dV- (A.6) 
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(Note that though the overlap of the ground states and as in Eq. (A. 3) is neghgible, 
the couphng coefficient K is not, due to the local maximum of the potential Kxt(p) in 
the overlap region, see Fig. 1 of section 1). Finally, assuming that Qu <^ (i-e., attractive 
atomic interaction in the w-condensate) and setting 



in Eqs. (A.4)-(A.5) one arrives at the dimensionless system (l)-(2) with a = —gv/du and 
/Xo = i^v — Su)/K, where the space and time variables of the system (l)-(2) are given as 
follows X = {\/2mK /h)x and t = {K/h)i. 

Recalling the well-known expression for the coupling coefficient gj^^^ = 4:7rh'^a'f-'> /m, where 
a(") is the scattering length in the -u-condensate, and using the definition of (?„ (A. 6) in the 
relation (A.7) we conclude that the actual (i.e., physical) numbers of atoms are related 
to the numbers of particles in the model (l)-(2) as follows 



where we have introduced the characteristic overlap distance by do = h/ \/2mK. 

In the above derivation it was assumed that the transverse profile of the order parameter 
is determined by the linear part of the r.h.s. in equation (A.l). Such approximation is 
valid for the weak atomic interaction: 

max(|gi"^|,|g^^^|)Ar nu;± 

where d± = {h/mw±Y^'^ and d\\ are the transverse and longitudinal sizes of the condensate, 
with the transverse size being given by the oscillator length, is the (physical) total 
number of atoms, and uj±_ is the external trap frequency (in the transverse dimensions) 
without the separation barrier. The r.h.s. of formula (A. 9) is an estimate on the energy 
of the hnear transverse terms in the GP equation (A.l). Substitution of the expression 

for the coupling coefficients g'o"'''^ — 47rfi^a^'"'"V"^ iii the condition (A. 9) leads to a more 
convenient equivalent form: 

87rmax(|a(")|, \a^^^\)^ < 1. (A.IO) 



To make an estimate we use the usual values: ~ 10 ^m and d\\ ~ 10 ^m (the latter 
corresponds to the asymmetry parameter 7 = d^/d\\^ ~ 0.01 for d±^ ~ 10~^m, usual for 
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experiments on the cigar-shaped condensates). Thus N <^ 10^ thereby allowing for up to 
a hundred thousands of atoms in the condensates. 

Due to the scaling d\\ ~ N, where N is the total number of particles in the model (1)- 
(2) (compare the widths of the solitons from figs. 3 and 4-5), the condition (A. 10) is 

satisfied for all values of the parameters fio and a in the system (l)-(2) and all values of 
the chemical potential fi if it is satisfied just at one point (a, /iq, /u). Therefore, though the 
total numbers of atoms do vary by an order of magnitude, nevertheless, the model system 
(l)-(2), in fact, applies uniformly in the whole parameter space. 

A very rough estimate on the numbers of trapped atoms can be easily carried out. Note 
that the coupling coefficient gu depends on the geometry of the trap through the multiplier 



(dP) 



2 



Here dj_' is the characteristic transverse size of the li-condensate. To make an estimate 
on the numbers of trapped atoms we assume: \a^^^\ ~ lO-'^m, dj*^ ~ = lO-^m, and 
do ~ O.lrf^"* = 10~^m, what leads to the transformation coefficients in formula (A. 8) 
of the order lO'"^. Thus this estimate shows that all the stable solitons we discuss in the 
present paper correspond to numbers of trapped atoms lying below the estimated bound 
of 10^. 



B Unstable mode and the related perturbation of the soliton solution 



Here we provide some details on the unstable eigenfunction of the hnear stability problem 
and the related initial perturbation. Looking for eigenmodes of the linear system (22) we 
write 







(x \ 












Xy 






[yJ 



since we need a real solution. The linear problem (22) then reduces to an eigenvalue 
problem for A = —iD,, which can be written as a system of two matrix equations: 







(yA 












= Ao 




, -in 























(B.2) 
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from which one derives equation (24). Note that is real (as the eigenfunction of 

the eigenvalue problem (24)). Then from the second equation in (B.2) we conclude that 
the eigenfunction (X„, Y^, Xy, Y^)'^ corresponding to the real eigenvalue (i.e. imaginary 
eigenfrequency fl) is real. Taking this into account and using equations (B.1)-(B.2) it is 
straightforward to deduce that the initial perturbation {ui,vi)^ of the soliton solution 
which is based on the unstable eigenfunction {X^, Y^, X^, Y^)'^ is given as 



where a is arbitrary real constant. 



C Negative eigenvalues of the operator Ai for a > 



To prove that Ai has only one nodeless eigenfunction (and, hence, only one negative 
eigenvalue) for a > 0, we use the fact that the in-phase soliton solutions correspond to 
the /ci-branch of the dispersion relation (11). The eigenfunction {Xi, X2Y of the operator 
Ai corresponding to an eigenvalue A satisfies the following system of equations, 



"^'^^ + (a^ + A + 3t/2) + - 0, (C.l) 
+ (/X + A - /xo - 2>aV^) X2 + Xi = 0, (C.2) 



dx^ 
dx^ 



which follow from the definition of Ai. We will prove that, among two possible combina- 
tions of nodeless eigenfunctions, X1X2 > and X1X2 < 0, the latter one is not possible 
for the in-phase solitons when a > 0. 

The multiplication of Eq. (C.2) by dX2/dx and integration from x = to x = 00 give 
(/i - /io + )^)X^{0) = Jdx l^-6aV''X2^ + 2^i^| • (CS) 



Here we have used that the eigenfunction has no nodes, hence dX2/dx = at x = 0. 
Taking into account that the condition > (see equation (11)) implies /i — fio < 0, 

and that X2(dX2/d,T) < for a; > 0, we conclude that Eq. (C.3) cannot be satisfied for 
a negative eigenvalue A if Xi(dX2/da;) > for x > 0, i.e., the operator Ai cannot have 
an eigenfunction corresponding to a negative eigenvalue and, simultaneously, having the 
property X1X2 < 0. 
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Distance 



Fig. 1. A scematic transverse shape of the trapping potential. The two condensates are designated 
by the symbols u and v. The difference of the chemical potentials between them is /xq- 
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Fig. 2. The numbers of atoms in the two condensates vs. the chemical potential fi for a = 
and HQ = 1, 0, —1, —5. The total number of atoms is given by solid curves, the number of 
atoms in the it-condensate by dotted and the number of atoms in the t;-condensate by dashed 
curves. Here the ?;-condensate contains non-interacting atoms, while the atomic interactions in 
the "u-condensate are attractive. 
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Fig. 3. The soliton solutions with a fixed total number of trapped atoms which correspond to 
the three values of chemical potential marked by stars in Fig. 2(c). Solid and dashed curves 
correspond to the u- and v- components of the solitons, respectively. 
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Fig. 4. The solitons with nearly all the atoms collected in the condensate with weak repulsion 
(a = 0.001). The left panel shows a zoomed-in part of the curves N = N{iJ,), Nu = N^in) 
and Ny = Ny{fi) about the local maximum. The soliton solutions corresponding to the stars 
on the curve A'" = N{ii) are given in the right panel (here the solid curves correspond to the 
-y-component, while dashed to the ti-component) . 
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Fig. 5. Typical curves for the number of atoms vs. chemical potential (left), and an example of 
the unstable soliton solution (right) featuring a large share of atoms in the repulsive condensate. 
Here a = 0.1 and /xq = —2. The soliton corresponds to the star on the curve N = Ar(/x) (the 
solid curve corresponds to the n-component and dashed to the v-component) . 
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Fig. 6. The number of atoms in the condensates vs. the chemical potential (top), and examples 
of the solitons (bottom) in the case when the soliton bifurcation from the zero solution is 
accompanied by a jump in the number of atoms (here a = 1 and /xq = 0). 
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Fig. 7. The bifurcation diagram for a = —0.1 and /io = \J—a, — a ~ —2.846. The top 

panel shows the total numbers of atoms in the condensates for the two branches of the soliton 
solutions. The bottom plot shows the numbers of atoms in each of the two condensates. In the 
figure, the thick curves correspond to the family of the sech-type solitons. 
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Fig. 8. The same as in Fig. 7 for a = —0.8 and /xq = —0.2, except that here /xq 7^ Mo(a)- 
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Fig. 9. The panels (a) and (b) show the bifurcation diagrams for a = —1 and /xq = 0, and the 
panels (c) and (d) show the diagrams for /xq = —0.01. As well as in Figs. 7 and 8, the top 

and bottom plots display, respectively, the total number of atoms, and the numbers of atoms 
in each condensate. The thick curves correspond to the sech-type branch of solitons and its 
deformations. 




Fig. 10. The unstable soliton and the real part of the single eigenmode of small perturbations 
which gives rise to instability. 
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Fig. 11. Evolution of the unstable perturbed soliton. In the top panel, the general direction 
of the evolution is indicated relative to the curve N = N{iJ,). The cases 1 and 2 correspond, 
respectively, to the addition and subtraction of the unstable mode (the one shown in Fig. 10), i.e., 
to the small perturbation which equals the unstable eigenmode multiplied by a small positive or 
negative amplitude. The bottom panel: the corresponding evolution of the numbers of particles 
in the two condensates. 
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Fig. 12. The decay of an unstable soliton into radiation in the case when dN/d^ ^ oo as 
A* ~^ Mmax- This evolution corresponds, for instance, to the stability curve shown in the left 
panel of Fig. 5. 
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